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Abstract 

Fluid-structure-interaction (FSI) is a very important 
consideration in the modeling of fluid flow over and 
through contemporary aircraft and aerospace vehicles. The 
numerical code used in this research utilizes a Lagrangian 
particle-based solver in conjunction with an Eulerian mesh- 
based solver (known as the Material Point Method (MPM)) 
to gain the advantages that both solvers offer to reduce 
computational overhead. Heat transfer between a fluid and a 
solid is an important FSI consideration, especially for high 
velocity (supersonic and hypersonic) aircraft and aerospace 
vehicles. A least squares gradient calculating technique was 
utilized combined with Fourier's Law of Conduction to 
incorporate heat transfer into the MPM algorithm. 
Application of the heat transfer algorithm to hypersonic 
single ramp and double cone flows was investigated. The 
MPM results were compared to experimental and numerical 
data from independent sources and observed trends agreed 
well with the results from these sources. 
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Nomenclature 

«o > #1 > #2 ' a 3 ' a A ' a 5 = me coefficients to be solved 
for in the least squares equations 
a = acceleration vector 
b = body forces vector 

D = particle diameter for Knudsen number 
calculations 

e = internal energy of material point 

e res = the residual error of the assumed solution fit 

S = strain 

S = strain rate 



f ex = external force vector 
f int = internal force vector 

F = the evaluated value for the assumed least squares 
solution fit 

Ggp = gradient of the shape function 

Y = ratio of specific heats 
I = second order unit tensor 
Kn = Knudsen number 
k = thermal conductivity 
kfr = Boltzmann constant 

L = characteristic length 
m = mass 

ju = absolute viscosity 

p = density of material point 

p = momentum 

p = pressure 

q = heat flux 

(J = Cauchy stress tensor 

S gp = shape function for transferring values between 

material points and grid 

S r = sum of squares of the residuals 

T = temperature of material point 

T s = temperature at solid surface (wall temperature) 

T x = temperature of ambient flow 
t = time variable 
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u c = calculated solution value 
u e = experimental solution value 

V = velocity vector 
X = position vector 

X = x-position variable used in least squares algorithm 

y = y-position variable used in least squares 
algorithm 

Introduction 

The Material Point Method (MPM) has drawn much 
interest in recent years for modeling solid and fluid 
dynamics problems due to the innate ability to 
monitor individual particle positions without loss of 
accuracy or difficulty in coding. Within MPM, the 
momentum equation is solved on an Eulerian mesh; 
velocity and position information is transferred back 
to the individual particles. The MPM code used in this 
research (referred to hereafter as the Aeroelastic 
Material Point Method (AEMPM)) is uniquely able to 
perform these calculations for both solid and fluid 
phases within a single numerical algorithm, which 
allows fluid-structure interaction in the momentum 
(force) equation. However, AEMPM does not 
currently include capability to transfer heat (energy) 
between the solid and fluid. The motivation for the 
current research is the development of an integral heat 
transfer algorithm within a fluid and solid particle 
based numerical solver. This would allow both 
momentum and heat transfer phenomena to be solved 
in a single numerical code without the need for 
numerical coupling between multiple solvers. 

This paper discusses a method to incorporate heat 
transfer into a MPM code and the results generated. A 
least squares gradient technique is used to calculate 
temperature gradients and a Fourier's Law of 
Conduction heat transfer algorithm utilizes this least 
squares result to calculate heat transfer between the 
fluid and solid. Aerodynamic heating by viscous shear 
and radiation heat transfer are not modeled. The 
hypersonic flows considered in this paper are laminar 
and so turbulence effects are not modeled. Verification 
studies of the algorithm with a 35° ramp and a 25°/55° 
double cone in a hypersonic flow of nitrogen are 
performed with numerical and experimental data. 

The MPM code utilized in this research is a Navier- 
Stokes based continuum solver. This type of code was 
selected in lieu of a Direct Simulation Monte Carlo 
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(DSMC) particle based code since the flows analyzed 
are in the continuum regime. Additionally, it is very 
impractical to use a DSMC solver on denser flows 
(Knudsen number less than 0.01) due to the large 
number of particles necessary to resolve the flow 
physics. 

State-of-the-art Review 

Numerical modeling of heat transfer is significant 
within a code designed for hypersonic vehicle analysis. 
Hollis et al. found the peak heating load for 
hypersonic flight of the X-33 aerospace vehicle to be 
2 

around 100,000 W I m . With such a large heat load 
on the vehicle, accurately modeling the heating effects 
during flight to determine material capabilities and 
thermo-mechanical fatigue effects is vital. 

In order to integrate heat transfer into a computational 
architecture, temperature values or gradients must be 
known. A least squares algorithm was integrated 
during the preliminary development of the AEMPM 
code (see multiple papers by Hu et al.) to calculate 
property gradients of velocity, density, and energy 
based upon techniques discussed in the paper by Koh 
and Tsai. Koh outlined a least squares curve-fitting 
technique and further defined the setup and numerical 
implementation of surface boundary conditions into a 
Lagrangian algorithm. 

During the development of a numerical algorithm, it is 
necessary to obtain (or generate numerically and 
experimentally) comparable cases that can be used to 
verify and validate the algorithm. For hypersonic 
flows, the single ramp and double cone are commonly 
used baseline cases in the literature (see Harvey et al., 
Gnoffo, Candler et al., Nompelis, and Moss et al.). 
Coleman and Stollery performed a study on a single 
ramp at hypersonic flow velocities and found heat 
transfer rates on the order of 500,000 to 1,000,000 

W I m . Moss et al. studied hypersonic heat transfer 
of both the single ramp and the double cone. The 
double cone was very similar in geometry to a single 
ramp (from a two-dimensional perspective) with the 
addition of a second ramp downstream and heat 
transfer rates were of similar orders of magnitude as 
the single ramp. 

Moss et al. utilized a Direct Simulation Monte Carlo 
(DSMC) for Mach 25 flow over a control surface (i.e., 
an elevator or aileron). The numerical simulations 
were for a 35 degree compression ramp at a low 
density flow condition, the results of which were 
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compared qualitatively with oil flow pictures. The 
DSMC numerical code, a general 2-dimensional 
(known as G2) or general 3-dimensional (known as F3) 
axisymmetric code, is a particle based method in 
which particles are populated in discrete cells. 
Molecular collisions were simulated with the variable 
hard sphere molecular model and energy exchange 
between kinetic and internal modes was controlled by 
the Larsen-Borgnakke statistical model. The surface of 
the single ramp was defined with a specified constant 
temperature and full thermal accommodation and 
diffuse reflection were assumed for the fluid-structure 
interactions. The single ramp results were found to be 
sensitive to grid resolution, however, these results 
showed good qualitative agreement with the oil flow 
pictures. 

Moss et al. utilized a revised version of the G2 code, 
known as the DS2V code, to analyze several double 
cone cases. The DS2V code models both time accurate 
unsteady flow and time-averaged steady flow 
simulations as well as molecular collisions and energy 
exchange in the same manner as the G2 code. Double 
cone cases were used for validation due to the 
availability of large amounts of experimental and 
numerical data. The double cone cases also exercised a 
numerical code's compression modeling capability 
since the double cone geometry introduced shock 
wave interactions at supersonic flow velocities. The 
results in Moss et al. found excellent agreement 
between the DS2V code and other numerical codes, 
however, an error of up to 15% was noted for heat 
transfer results when compared with experimental 
data. There was a single experimental case where the 
DS2V results were erroneous by up to 30%, but the 
experimental results were believed to be incorrect. 

Chanetz et al. analyzed a 25765° double cone in Mach 
10 flow conditions using two Navier-Stokes based 
continuum numerical solvers and two DSMC solvers 
(including the G2 code). These computations were 
compared to experimental data obtained in the R5Ch 
wind tunnel in which the mean free path was equal to 

5*10 ^meters. Chanetz et al. concluded that this 
mean free path justified the use of the Navier-Stokes 
equations as the experiments were performed in the 
continuum regime. The computational results between 
the Navier-Stokes continuum and DSMC solvers 
agreed and trended with the experimental data except 
for a 33% deviation at the maximum heating rate 
where the shock wave and boundary layer met 
(computational solvers predicted 33% low). The 



computational solvers' prediction of the separation 
location disagreed with experimental data by about 
17%. 

Harvey reviewed the state-of-the-art for the hollow 
cylinder/flare body (similar to an engine inlet) and the 
25755° hypersonic double cone cases for both Navier- 
Stokes and DSMC solvers. These cases were in the 
realm of Navier-Stokes solutions due to their relatively 
high freestream densities, but close enough to the 
rarefied realm that DSMC solvers could practically be 
applied, although, the computational demands were 
extremely high. The DSMC solvers required very fine 
resolution to resolve the flow and underpredicted 
separation lengths even with cell adaptation and 
collision modeling improvements. The Navier-Stokes 
solvers tended to exhibit spurious perturbations in the 
flow due to the flow's steep property gradients. 
Harvey concluded that the flows studied were 
insensitive to vibrational temperature fluctuations 
suggesting that the analyzed flows were vibrationally 
frozen. 

Material Point Method Algorithm 

The fundamental governing equations in the Material 
Point Method are the conservation equations for mass, 
momentum, and energy (equations 1, 2, and 3, 
respectively; see Wikipedia reference for a brief 
description of the Material Point Method). In AEMPM, 
the time derivatives were discretized using the 
backward Euler method and the spatial discretization 
was performed in a manner similar to typical Finite 
Element Methods using shape functions and the weak 
forms of the momentum and energy equations (see 
Bardenhagen and Kober). 



dp 

dt 



+ Vp-v = 



pa = Vcr + pb 



pe = a : s 



(1) 



(2) 
(3) 



The strain rate in equation 3 is defined in equation 4. 
For the present investigation, the small deformation 
assumption was valid and the strain rate relation in 
equation 4 was used in the calculations. 



a = 



1 



2 L 



(Vv) + (Vv) 7 



(4) 



These equations were discretized for material points, 
the grid nodes on the Eulerian mesh, or both. The 
index convention for the discretization equations 
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followed that in Hu et al. and were similar to York II et 
al., in which grid variables have a subscript g and 
material point variables have a subscript p . The 
remaining equations in this paper were derived from 
these two references. The momentum equation 
(equation 2) was discretized on the grid as shown in 
equation 5. 



m g a g =f 



int f ext 



(5) 



Equation 5 was modified to a form based on the 
relation between acceleration and momentum in 
equation 6. Next, updated values of momentum at the 
next timestep were calculated as shown in equation 7. 



— — = ma 

dt 88 



(t + At) = p g (t) + At{f g * + fl x \ 



(6) 



(7) 



material point stresses were calculated as shown in 
equation 13. 

Ep {t + At)= Ep (t)+^Y J l G 8P V 8 ( f + Af ) + ( G 8P V 8 ( f + At f 1 



(12) 
(13) 



Since the external force vector was an input within the 
AEMPM code, only updated values of the internal 
force vector must be calculated for use in the 
momentum equation. The internal force vector was 
calculated using the material point stresses in equation 
14. Then, the fluid material point density was updated 
for continuity (equation 1) using equation 15. Note 
that the value inside the trace (the summation of the 
main diagonal of a matrix) of the denominator is the 
strain increment for the fluid. 



The above discretization equations required values to 
be located on the grid; however, the solution domain 
was initialized with material points, so shape 
functions (see Bardenhagen and Kober for a 
description of the shape functions used) must be used 
to transfer information between the grid and material 
points. Values of mass and momentum were 
calculated using the shape functions in equations 8 
and 9. 



m 



=Z 

P 



m p S gp 



m g y g 



= 1 

P 



m p^ p^ gp 



(8) 



(9) 



The material point velocity and position were then 
solved using the grid values of acceleration and 
momentum as shown in equations 10 and 11. 



v p (t + At) = v p (t)+Y J S gp a g At 

g 



(10) 



• int 



m, 



8P P P P 



(14) 



p p {t + At) = - i Pp pfl (15) 

/ 1 + tr[s p (t + At) -£ p {t)J 

The material point energy was updated using the new 
values of density and stress in equation 16 and the 
material point pressure was calculated with updated 
values of density and energy using the equation of 
state as shown in equation 17. The material point 
temperature was calculated by dividing the material 
point energy by the specific heat of the material for 
both fluid and solid phases. 



eJt + At) = eJt) + 



p P =ir- 1 )p P e 1 



a p (t + At\s p (t + At)- s p {t)) 



p p {t + At) 



(16) 
(17) 



x p {t + At) = x p {t) + Y,S gp 1 ^ L At + ^ S gP a g( At ) 2 

g m g g 

(11) 

The material points have values of updated positions 
and velocities which were used for calculations of the 
updated strain values (for the solid material) using 
equation 12. The term summed with the strain value at 
time t is known as the strain increment. The fluid 



Heat Transfer Algorithm 

This paper focuses on the fluid-structure interface and 
the transfer of thermal energy across this interface. 
The computational layout of the AEMPM algorithm 
used in this research is shown in Fig. 1. The grid 
remains fixed throughout the entire calculation and 
was used as a means of momentum transfer. The 
material points were initiated in the solution domain 
with mass and flow properties. Figure 1 shows the 
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separate types of material points present in the 
algorithm. The fluid and solid MP's represent the fluid 
and solid phases respectively, while the surface MP's 
represent the fluid-structure interface. These surface 
MP's, initiated in the algorithm with fluid properties, 
were used as the exchange "medium" for energy 
between the fluid and solid MP's. A "cutcell" is 
defined as a grid cell in which the fluid-structure 
interface crosses. For the hypersonic single ramp and 
double cone cases, the Dirichlet boundary condition of 
a constant wall temperature was enforced to allow 
comparison with experimental data and other 
numerical results that enforced this thermal boundary 
condition. The solid MP's were omitted from the 
analyses since a constant wall temperature boundary 
condition yielded no temperature change within the 
solid material. 



Non- 
CutCell 



CutCell 



Surface MP 




Solid MP 



Fluid-Solid 
Interface 

FIG. 1. AEMPM COMPUTATIONAL LAYOUT. 

In the present study, as a fluid flows over a ramp and 
double cone at high speeds, large pressure gradients 
(strong compression) are developed in the flow which 
result in strong shock waves, leading to large 
temperature changes of the fluid. The temperature 
distribution in the interface region is continuous (i.e., 
there is no temperature jump). However, the fluid 
particles in contact with the solid wall may not be the 
same temperature as the solid because there is a small 
distance between the fluid and solid particles where a 
temperature gradient exists. This temperature 
difference decreased as grid resolution increased and 
was less than 5% of the wall temperature for the 
highest grid resolution cases. Viscous shear also 
generates heat transfer, but this friction heating 
mechanism is not currently modeled in the numerical 
algorithm and so translational energy (fluid 
compression) is the only mode of fluid particle 
temperature change. 

Fourier's Law of Conduction for a Cartesian grid is 
shown in equation 18 (from Incropera et al.). The heat 



flux between two points can be calculated using the x- 
and y- derivatives of the temperature (calculated using 
the derivatives from the least squares solution in 
Appendix I), which is shown in equations 19 and 20. 
The ' a ' coefficients in equations 19 and 20 are the 
same coefficients as those seen in the least squares 
solution from Appendix I. 



= -k 



8T 8T 

— + — 
dx dy 



dx 

dT 

~dy~ 



a.\ + la^x + a^y 



= a-2+ 2a^y + a$x 



(18) 



(19) 



(20) 



The boundary conditions utilized in the MPM results 
in this paper were characteristic boundary conditions 
for the inlet, outlet, top, and bottom fluid domain 
boundaries (see Anderson Jr. for characteristic 
boundary condition description). The fluid-structure 
interface was defined as a Dirichlet boundary 
condition in which the temperature of the solid surface 
was fixed. 

Application of Heat Transfer to a Single 
Ramp in Hypersonic Flow 

The first verification case was a single ramp at a 35° 
angle. A single ramp is a simplified model of a 
hypersonic vehicle's control surface (such as an 
elevator) and exercises a numerical code's 
aerodynamic flow and compressibility simulation 
capabilities. A motivation to simulate these types of 
flows is driven by the fact that flight safety greatly 
depends on the effectiveness of control surfaces. 
Numerical data from Moss et al. has been compared to 
the results from the AEMPM code in this paper. The 
35° single ramp was setup as shown in Fig. 2 with the 
coordinate locations in meters (L = 0.0714 meters). The 
flow was nitrogen with constant specific heats and had 
the following initial conditions: =1521 mis , 

^ =1.401 *10" 4 kg/m 3 , 7^=9.06 K , 

= 0.376 Pa , Moo = 24.8 , and T s = 403.2 K . 

Note that the surface of interest is made up of 
coordinate locations A, B, and C (the heat transfer 
surface of a flight control surface deflected in the 
positive y-direction) and the results shown reflected 
only this surface. 
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A 


B 








E 




L 





Point 


Coordinate 


A 


(0,0) 


B 


(0.07140,0) 


c 


(0.12988.0.04695) 


D 


(0.12988.-0,006) 


E 


(0.02239,-0.006) 



FIG. 2. SINGLE RAMP CONFIGURATION AND SETUP 
COORDINATES. 

The Knudsen number was calculated using equation 
21. For the single ramp flow condition, the Knudsen 
number was 0.007 and therefore, the continuum 
assumption (Navier-Stokes equations with the no-slip 
boundary conditions) is valid (Panton). For reference, 
Panton classified the flows as follows (Panton, page 
680): 

0<i&z<0.01 "continuum" flow (Navier-Stokes 
equations with no-slip boundary conditions) 

0.01 < Kn < 0.1 slip flow (Navier-Stokes equations 
with slip boundary conditions) 

Q.\< Kn<?> transitional flow (Navier-Stokes 
equations are not valid in this regime) 

3 < Kn < oo free molecular flow (molecules 
interact only with walls) 



Kn 



4lnD 2 pL 



(21) 



The results of the MPM code and the verification case 
are shown in Fig. 3 (zoomed in at the ramp transition 
from x=0.05 m to x=0.13 m). The verification case is 
another independently developed numerical code 
(known as the G2 code) which is discussed in more 
detail in Moss et al. The number following the "MPM 
Grid=" in the legend is the relative resolution of the 
grid used in the computational domain. Grid 
resolution and the number of material points per grid 
cell are both variables within the AEMPM code. 
During the studies discussed in this paper, it was 
found that the optimal material point density was 4 
material points per grid cell (2x2 pattern in x- and y- 
directions). Additional material points per grid cell 
beyond this caused instability and larger perturbations 
in the numerical algorithms. The results discussed in 
this paper utilized four material points per grid cell 
and varying grid resolutions. 



The differences between the G2 code and MPM results 
are attributable to the differing code formulations. The 
G2 code, a Direct Simulation Monte Carlo numerical 
code, includes energy exchange between translational, 
rotational, and vibrational modes. Moreover, diffuse 
reflection is assumed for the fluid-structure interface. 
The AEMPM code only incorporates the translational 
energy exchange mode and is a Navier-Stokes 
continuum based numerical solver which excludes 
diffuse reflection. 




FIG. 3. SINGLE RAMP HEAT TRANSFER DISTRIBUTION. 
2 

The L relative error norm was calculated for the 
MPM data using equation 22 and is shown in Table 1. 
Increasing the grid resolution by about 50% (from 
MPM Grid=270 to MPM Grid=400) reduced the 
relative error norm by about 60%. Note the large 
reduction in relative error for the peak heat transfer 
rate (near x=0.12 m in Fig. 3). 



L error = 



(22) 



TABLE 1. MPM HEAT DISTRIBUTION COMPARISON TO 
NUMERICAL RESULTS. 





2 


Relative Error in 




L Relative 


Peak Heat 




Error Norm 


Transfer Rate 


MPM Grid=270 


0.265 


36% 


MPM Grid=400 


0.101 


4% 



Application of Heat Transfer to Double Cone 
in Hypersonic Flow 

The double cone is a popular case for studying 
hypersonic flows (a few different configurations and 
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cases were discussed in Harvey et al.) as it 
incorporates complex shock wave interaction between 
two shocks and the associated heat transfer at the 
fluid-structure interface. For this investigation, the 
double cone was setup as shown in Fig. 4 (L = 0.092072 
meters). The coordinate locations are in meters and the 
points below the centerline are the same x- values with 
a negative y-value. The first cone has a half angle of 
25° and the second cone has a half angle of 55°. 




Point 


Coordinate 


A 


(0,0) 


B 


(0,092072,0,042934) 


C 


(0.153*66,0.130900) 


O 


(O.195OO0.O.1309O0) 



FIG. 4. DOUBLE CONE CONFIGURATION AND SETUP 
COORDINATES. 

The case studied had the same initial conditions as 
Moss et al. (nitrogen at = 2072.6 mis , 



Poo 

Pr, 



= 1.757*10" 
2.23 Pa , M c 



kg I m~ 



42.61 K , 



15.6 , and T. = 297.2 K 



and constant specific heats. The Knudsen number for 
the Double Cone flow conditions was 0.004 and so the 
continuum assumption is valid. The AEMPM results 
in this paper were compared to the experimental 
values from Moss et al.'s paper for increasing degrees 
of grid resolution in Fig. 5. 



140 - 



,-120 
E 

§100 



80 - 
60 - 
40 P 

20 

c 






Exper. Data (Moss 2005) 




MPM Grid=20D 




MPM Grid=300 




MPM Grid=50D 




0.05 0.1 

x(m) 

FIG. 5. DOUBLE CONE HEAT TRANSFER DISTRIBUTION. 



The results exhibited similar trending as the 
experimental results from Moss et al.'s paper. The 
double cone MPM results did not exhibit the "dip" 
trend in the separation region (the area between x=0.07 
m and x=0.1 m) because the AEMPM code does not 
resolve recirculation or flow separation effects. 
However, note that the peak heat transfer 'spike' seen 
at x=0.1 m is resolved within about 6% in the most 
resolute case. Furthermore, note the oscillations seen 
in the results which correlate with the observations by 
Harvey that Navier-Stokes solvers tend to exhibit 
spurious perturbations due to large property gradients 

2 

The L relative error norms and relative errors in 
peak heat transfer rate calculated for the MPM data 
are shown in Table 2. Note that increasing the grid 
resolution from 200 to 500 (approximately six times 
the MP's in the computational domain) reduced the 
relative error norm by about 65% and significantly 
improved the relative error in the peak heat transfer 
(near x=0.1 m). 

TABLE 2. MPM HEAT DISTRIBUTION COMPARISON TO 
EXPERIMENTAL RESULTS. 







2 

L Relative 
Error Norm 


Relative Error in 
Peak Heat 
Transfer Rate 


MPM Grid= 


=200 


0.584 


63% 


MPM Grid= 


=300 


0.443 


48% 


MPM Grid= 


=500 


0.207 


6% 



Conclusions 

The incorporation of heat transfer within a hypersonic 
aerospace vehicle numerical model is vital to calculate 
the effects of the temperature gradients and 
compression developed within shock waves and 
pressure gradients at the fluid-structure interface. The 
least squares gradient method described in this paper 
was used to calculate the temperature at the fluid- 
structure interface and heat transfer at this interface. 
Using a conductive model (Fourier's Law of 
Conduction) for fluid-structure interface heat transfer 
resulted in agreeable trends for the single ramp case 
and a preliminary study in the double cone case. 

ACKNOWLEDGMENTS 

Funding for this research and the base numerical 
AEMPM code were provided by Advanced Dynamics 
Incorporated, Lexington, Kentucky. 



112 



Frontiers in Aerospace Engineering Volume 2 Issue 2, May 2013 



www.fae-journal.org 



Appendix I 

Least Squares Gradient Technique 

The least squares gradient method is a matrix algebra 
solution technique used to calculate gradients of a 
function value. The AEMPM code uses a least squares 
gradient technique to calculate values of velocity, 
density, and energy with a Dirichlet boundary 
condition (known value on the boundary) for velocity 
and energy as well as a Neumann boundary condition 
for density. This paper studied the implementation of 
the least squares gradient technique into the AEMPM 
code to calculate temperature gradients using the 
temperature information of neighboring material 
points. The basic steps used to implement the least 
squares method for temperature gradients have been 
outlined below. A second order polynomial was 
utilized due to its ease of implementation and yielded 
more accurate results over a first order polynomial. 
Other functions such as cosine series were investigated, 
but stability issues were observed and so a second 
order polynomial was utilized for the implementations 
discussed in this paper. 

Assume a solution of the form: 

2 2 

F = ciq + a\x + a^y + a^x + a^y + a^xy + e res 

(AI.l) 

Solving this equation for the sum of the squares of the 
residuals e res by means of the least squares method 
yields (Chapra and Canale): 

n 

S r =XX F i -a -a l x j -a 2 yj - a sXj 2 - a 4 y t 2 - a 5 x iyi ) 2 
(=1 

(AI.2) 

To solve for the unknowns, the derivatives of this 
equation are taken with respect to each unknown 
coefficient yielding the following set of equations: 



dS r 
da. 



= - 2 Z( F ' " a o ~ a \ x i ~ a iyi - a 3 x i - a 4yt - a 5 x iy0 



f=l 



dS n 

—^ = -2^x i (F i -a Q -a x x t -a 2 y i -a 3 x i 2 -a 4 y 2 - a 5 x ; ;y,) 



dS r 
da 



= - 2 ^yi( F i-ao-aix i -a 2 yi-a 3 x i -a 4 yj -a 5 x iyi ) 



2 i=\ 



(AI.3) 



dS r 
da 



9S r -v> 2,,, 2 2s 

- — = -22 J y i {F i -a -a l x i -a 2 yi-a 2l x i -a 4 y t -a^yO 
3a 4 , =1 



dS r 
da 



- r ^Hyii F i-a< A -a l x i -a 2 y i -a- i x i -a 4 y t -a 5 x iyi ) 



5 <=1 

Finally, setting these equations equal to zero to find 
the minimum of the sum of the squares and re- 
arranging them yields the following system of 
equations (where all summations are from i=l to n, 
with n being the number of material points included in 
the solution): 

(n)a +(2>,)ai +(X> - ;)«2 +[^x i 2 )a 3 +(X>'< 2 ) a 4 + (Z ^ )"5 =Z 



(Z *> )"o + (z x i 1 ) a i + (Z x i y- )°2 + (Z *> 3 ) a 3 + (Z x i y< 2 )°4 + fe x ; 2 y> K = X x < 

(Z y< k> + (Z x i y< K + (Z y> 2 ) a i + (Z *• 2 y< ) a 3 + (Z y> 3 ) fl 4 + [Z ^ >'< 2 ) a 5 = Z 

(AI.4) 

E^' 2 ^ +fZ x i 3 ] fl i + vL x i 2 yi\ , 2 +fZ ;c i 4 l 3 + [ S*i 2 yi W +(Z x i 3 >f) 2 5 =Z x i 2 



2 2^ 



fZ>'i 4 }'4+E :c i}'«- 3 } j 5=Z3 



= ^ 2 X x i ( F i ~ a o - a i x i ~ a 2yi ~ a i x i ~ a 4yi -a=,x iyi ) 



3 i=l 



(X-^iK + C Jc i 2 3'l7 ! l +(E x >'>'| 2 )°2 +[Z x >' 3 >'i) a 3 +(Z x i>'i 3 ) a 4 +(Z x r' 2 >'i 2 ) :, 5 =Z x i>'' 

Note that this system of equations is a symmetric 
matrix. Therefore, in the solution algorithm only a 
vector of the upper diagonal matrix needs to be 
generated. This vector can be generated in a loop very 
quickly by adding each material point's respective 
value for the summation in the upper diagonal being 

solved for (i.e., X,; ) or ^ J). The right hand 

side of equation AI.4 is solved in the same manner 
with summations in a loop. Finally, Cholesky 
Decomposition (see equations AI.5 and AI.6) is used to 
solve the system of equations and the original value of 
F from equation AI.l above can be obtained. This F 
is the least squares curve fitted value for a variable at 
the material point (i.e., temperature in the case 
discussed in this paper). 

Cholesky Decomposition 

Cholesky decomposition is a technique for solving 

symmetric matrices ([A] = [A]^j which can be 

T 

decomposed as [A]=[L][L] . To obtain the 
individual elements within the [L] matrix, equations 
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AI.5 and AI.6 are used (where equation AI.5 is for i 
indices from 1 to k — 1 ) from Chapra and Canale. 

i-l 

a ki ~ ^Tjhjhj 

l ki = " (AL5) 

Hi 



hk= ] a kk-Yu l lj ( AL6 ) 

V 7=1 
Appendix II 

AEMPMHeat Transfer Algorithm Pseudo Code 

The integration of the Fourier's Law of Conduction 
heat transfer algorithm for the fluid-structure interface 
into the AEMPM code followed the subroutine steps 
seen in the pseudo-code below: 

If cell is a cutcell 

For each MP in cutcell 

Loop through a 10x10 grid stencil centered 
around cutcell 

(loop includes all 100 grid cells in this step) 

Calculate relative position of fluid 
MP's and surface MP's in cutcell 

Calculate least squares coeff. 
with surface and fluid MP's 

Solve for unknowns using Cholesky 
Decomposition 

Calculate heat flux at surface MP using 
equation 18 
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